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ABSTRACT 

We present a new publicly available tool (DustPol) aimed to model the polarised thermal dust emission. The module DustPol, 
which is publicly available, is part of the ARTIST (Adaptable Radiative Transfer Innovations for Submillimetre Telescopes) package, 
which also offers tools for modelling the polarisation of line emission together with a model library and a Python-based user interface. 
DustPol can easily manage analytical as well as pre-gridded models to generate synthetic maps of the Stokes /, 2, and U parameters. 
These maps are stored in FITS format which is straightforwardly read by the data reduction software used, e.g., by the Atacama Large 
Millimeter Array (ALMA). This turns DustPol into a powerful engine for the prediction of the expected polarisation features of a 
source observed with ALMA or the Planck satellite as well as for the interpretation of existing submillimetre observations obtained 
with other telescopes. DustPol allows the parameterisation of the maximum degree of polarisation and we find that, in a prestellar 
core, if there is depolarisation, this eff'ect should happen at densities of 10^ cm~^ or larger. We compare a model generated by DustPol 
with the observational polarisation data of the low-mass Class object NGC 1333 IRAS 4A, finding that the total and the polarised 
emission are consistent. 

Key words, radiative transfer - methods: numerical - polarisation - submillimetre: general 



1. Introduction 

The importance of magnetic fields in regulating the physics of 
star formation is still a matter of debate. The main schools of 
thought invoke either the leakage of magnetic support through 
ambipolar diflTusion (Mestel & Spitzer 1956) or the supersonic 
turbulence which compresses the gas at the interfaces between 
converging flows (Elmegreen 119931) . While turbulent energy 
is required to sustain a cloud over its lifetime (Gammie & 
Ostriker I1996K ohmic dissipation time of magnetic fields is 
much larger than the collapse time, assuming a typical ionisa- 
tion fraction of ~ 10"^ - 10"^ and density of ~ 10^ cm"^ (Pinto 
& Galli I2Q08I) . Collapse time scale for a subcritical core is a 
few million to 10^ yr assuming the ambipolar diflTusion to be 
the dominant process (Mouschovias et al. 2006 ). Conversely, if 
collapse is controlled by turbulence, the time scale is shorter, of 
the order of 10^ yr (Mac Low & Klessen .2004j . Cores can also 



be forced into collapse because of external events such as col- 
lision between clouds as well as supernova shock waves with a 
time scale as short as 10^ yr (Boss et al. 2010). All these diflTer- 
ent mechanisms most likely take part in cloud collapse, but the 
role of magnetic fields cannot be ruled out from star formation 
theories, even when a core is supercritical (Basu il997i) . 

From an observational point of view, it is possible to gather 
information about the strength and the geometry of magnetic 
fields through measurements of the Zeeman splitting of hyper- 
fine molecular transitions as well as measurements of polarisa- 
tion of thermal dust emission in the (sub)millimetre wavelength 
range. Observations of the Zeeman eflTect provide information 
about the line of sight component of the magnetic field and 
have been carried out using atoms and molecules which show 
large magnetic moments (e.g. Crutcher et al. 1996, Falgarone et 
al. l2QQ8b . By observing maser emission it is possible to detect 
both circular and linear polarisation, making it possible to de- 
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rive the full 3D magnetic field (Vlemmings et al. 20111 Surcis et 
al. l2Qllall2Qllbb . On the other hand, polarisation of the thermal 
dust emission probes the plane of the sky component of mag- 
netic fields, providing information about the field geometry in 
molecular clouds and filaments associated with star formation 
regio ns (e.g . Novak et al. ll989 , Matthews et al. l2QQll Davidson 
et al. "2011 ). (Sub)millimetre wavelength interferometric obser- 
vations, first with the Berkeley-Illinois-Maryland Association 
array, BIMA (e.g. Girart et al. [19991 Lai et al. [20Q3]) . and nowa- 
days with the S ubmillimete r Arra y, SMA (Girart et al. i2009. 
Rao et al. I2QQ91 Tang et al. 120091) have proven useful to study 
the magnetic field properties at high angular resolution, probing 
scales of a few hundreds AU at which the dynamical collapse of 
the star forming cores is occurring. 

The Atacama Large Millimeter Array (ALMA) is the 
largest interferometer and the largest ground based project in 
(sub)millimetre astronomy. It provides a considerable enhance- 
ment in resolution, sensitivity and imaging with respect to pre- 
vious interferometers operating at these wavelengths. Since dust 
and a multitude of spectral lines can be observed at the same 
time, ALMA allows an in-depth investigation of the physical 
and chemical properties of e.g. embedded protostars, zoom- 
ing in to AU scales of nearby star- forming regions. ALMA 
will also open new possibilities for studying magnetic fields. In 
fact, due to the receiver capabilities, full polarisation calibra- 
tion and imaging will be the norm rather than the exception. 
For these reasons, providing ALMA users with self-consistent 
tools for polarisation modelling will be very important. The 
publicly available Adaptable Radiative Transfer Innovations for 
Submillimetre Telescope^ (ARTIST) package (J0rgensen et al. 
in prep.) provides a direct link between the theoretical predic- 
tions and the quantitative constraints from the submillimetre ob- 
servations, being able to model the full multi-dimensional struc- 
ture of, e.g., a low-mass protostar, including its envelope, disk, 
outflow, and magnetic field. 

The model suite is based on LIM^(L ine Modelling Engine, 
Brinch & Ho gerheij de 120 1 Ob , a non-LTE radiative transfer code 
using adaptive gridding in three dimensions that allows simula- 
tions of sources with arbitrary multi-dimensional structures en- 
suring a rapid convergence. In addition, ARTIST ofl'ers unique 
tools for modelling the polarisation of dust and line emission and 
an extensive library of theoretical (semi-analytical) models. All 
these packages are connected by a comprehensive Python-based 
graphical user interface with direct link to, e.g., ALMA data re- 
duction software (Common Astronomy Software Applications 
package, CAS^. 

In this paper, we describe the dust polarisation module of 
ARTIST: in Sect. |2] we briefly describe the polarisation of stellar 
radiation due to dust grains, in Sect. |3] we explain the features 
of the module for analysing dust polarisation emission, and in 
Sect. |4] we demonstrate the capabilities of the routines to handle 
analytical models, numerical simulations, and providing a com- 
parison to observations. Finally, in Sect. [5] we summarise our 
results and conclusions. 

2. Origin of dust polarisation 

Dust grain alignment provides a natural explanation for the 
origin of polarisation in the ISM. However, what mecha- 
nisms are responsible for grain alignment and what role 

^ http://youngstars. nbi -dk/arti St 

^ http : //www . nbi . dk/~br inch/lime . html 
^ |http : //casa . nrao . edu 



they play in diff'erent environments remain poorly understood 
(Lazarian 120031 12007b . The main explanations are supplied by 
paramagnetic (Davis & Greenstein 1951 ) and mechanical align- 
ment (Gold 1951 ), but other mechanisms have been proposed as 
well (Martin 1971 Purcefl 1979. Draine & Weingartner .l996j . 
Recent studies for example show that radiative torques, which 
are the results of the interaction of irregular grains with a flow 
of photons, can efl'ectively align the dust grains (Dolginov & 
Mvtrophanov lT976l Lazarian 120071) . 

Grains are most likely non-spherical since they are formed 
by stochastic growth processes leading to irregular shapes 
(Fogel & Leung 1998 ). Unpolarised starlight passing through 
an aligned-grain set will acquire polarisation due to the difl'er- 
ent grain extinction parallel and perpendicular to the alignment 
direction, since they absorb more light along their longer direc- 
tion. This depends on the fact that extinctions are proportional 
to the parallel and perpendicular grain cross section (see Lee & 
Draine 119851 and Efroimsky 2002 for a detailed discussion for 
spheroidal grains and a more general case, respectively). 

Conversely, in the submillimetre wavelength regime, the 
electric-field vector of the thermal continuum emission emerg- 
ing from rotating aligned grains has its maximum value in the 
plane containing the longer grain axis. The partially linearly po- 
larised vectors are parallel to the mean grain orientation, there- 
fore perpendicular to the magnetic field as projected on the plane 
of the sky. Consequently, from the observation of the polarised 
dust emission it is possible to model the morphology of mag- 
netic fields in the plane of the sky, the degree of polarisation 
depending on various contributions as the grain cross sections, 
the temperature, the composition, and the density. 

3. Dust polarisation tools 

With the dust polarisation module (DustPol) of ARTIST we aim 
to model the polarised dust continuum emission. The DustPol 
module is an extension to the LIME radiative transfer code. 
LIME is a non-LTE molecular excitation and line radiation trans- 
fer code that is mainly used to calculate line profiles in the far- 
infrared and submillimetre regimes. However, LIME can also 
ray-trace given density and temperature profiles and from that 
estimate the emerging continuum flux. It is this feature that we 
make use of in the DustPol module presented in this paper. At 
the relatively long wavelengths at which LIME operates, scat- 
tering is negligible and therefore the problem of making an im- 
age of the continuum model is reduced to integrating the source 
function 

h= B,{T^,,,)e-'^dT, (1) 
Jo 

along lines of sight through the model, ^^(rdust) is the Planck 
function evaluated at the dust temperature r^ust and dr^ = 
ka Pdust d^, where kx is the wavelength dependent dust opacity 
and Pdust is the dust mass density. The main benefit of using the 
LIME code for our application is its ability to automatically map 
a 3D model onto a computational grid. The grid points are placed 
at random in 3D space with probability of a point being placed at 
a given location weighted by a user- selected function, typically 
the underlying model dust density distribution. This results in a 
higher spatial grid point density in high dust density regions and 
vice versa for low density regions, which means that the expec- 
tation distance between neighbour points scales with the photon 
mean free path. 

Once the grid points have been placed, the Voronoi diagram 
of the distribution and the topological equivalent, the Delaunay 
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triangulation, is calculated (Ritzerveld & Icke 120061) . A Voronoi 
region associated with a certain grid point is defined as the set 
containing points that is located closer to the generating grid 
point than to any other grid point. The result is a tessellation 
of space consisting of random polyhedra. All physical proper- 
ties are taken to be constant over each polyhedra (or Voronoi 
cell). If the model resolution turns out insufficient in specific re- 
gions, more grid points can be added to these and the Voronoi 
diagram is recalculated. The grid is then regularised by moving 
each grid point slightly toward the centroid of the corresponding 
Voronoi cell several times in an iterative process. The result is a 
random grid consisting of approximately regular polyhedra (see 
Brinch & Hogerheijde 2010 for details). Ray-tracing is done by 
integrating straight lines of sight (rays) through the grid simply 
by summing up the constant contribution between the two inter- 
sections between a ray and the Voronoi facets for each Voronoi 
cell the ray passes through. 

For the dust polarisation calculations, the main extension to 
LIME is to include Stokes Q and U in the source function. We 
followed the procedure developed by Lee & Draine (11985) , then 
by W ardle & Konigl (1990), Fiege & Pudritz (200Q1), and Padoan 
et al. fSOOH) . As in Gongalves et al. ( 2005 ) and Frau et al. dSoTl ) . 
hereafter GGW05 and FGGll, respectively, our code accounts 
for the dependence of the emerging continuum radiation on dust 
temperature distribution (r^ust) supplied by the user (e.g., calcu- 
lated self-consistently through tools such as RADMC/RADMC- 
3E|^ Dullemond & Dominik 2004 ). According to this method, 
the polarisation degree p and the polarisation (or position) angle 
X, that is the direction of polarisation in the plane of the sky, are 
given by 



and 



1 n 
= 2 Jo 



max ^i(7^dust) ^ COS y - - dr^ , 



(8) 



(2) 



where a^max accounts for grain cross sections and alignment 
properties. For polarised grains, the eff'ective cross section is 
given by the diff'erence between the extinction cross section and 
a contribution coming from polarisation. For grains without po- 
larisation properties. Stokes / is basically given by Indeed, the 
grain properties may vary with density (Fiege & Pudritz 2000 ), 
but in the absence of detailed information, unless otherwise 
specified, we assume Ofmax to be constant and equal to 0.15 in 
the remainder of this paper, corresponding to a maximum polar- 
isation degree of about 15%, as suggested by GGW05. Higher 
values has not been observed and lower values appear unable to 
reproduce observations of the magnetic field of prestellar cores 
(FGGll). However, it is important to notice that the DustPol 
module allows the maximum polarisation degree to vary accord- 
ing to a user defined function (e.g. density or temperature depen- 
dent). 

After the ray-tracing, the code creates a file in FITS format 
which consists of three slices, one for each Stokes parameters. 
This FITS file can be read by the sim_observe-sim_analyze 
tasks in the CASA package to prefigure the expected polarised 
flux intensity map obtained by ALMA. 

4. Applications 

In this Section, a few applications of the DustPol routines are 
presented, showing their capabilities to manage input from ana- 
lytical models as well as numerical simulations. A comparison 
with SMA available observations is also provided. 



and 



tan 2;^ = 



(3) 



respectively, q, u, and / are called "reduced" Stokes parame- 
ters and are defined in terms of the standard Stokes parameters 
through a constant which depends on the cross section, a polar- 
isation reduction factor due to imperfect grain alignment and to 
the turbulent component of the magnetic field, and the Planck 
function. The "reduced" parameters are described by 



Q ~ I ^max 

Jo 



^^(^dust) e cos2j/r cos ydr^ 



U = I Qfn 

Jo 



^^(^dust) e ^mliff cos ydr^ , 



(4) 
(5) 



and 



4.1. Toy models 

It is instructive to start by considering geometrical models to 
check the eff'ective accuracy of the code. We examine sources 
with a purely toroidal (Btor), a dipole (Bdip), and a quadrupole 
(Bquad) magnetic field whose spherical components (Br, Bq, B(p) 
are 



Btor cx (0,0, r-^), 

f 2cos^ Vl - cos2 6 

Bdip oc 
and 



, 



, 3(3cos2£(-l) 3COS0VT - COS^ u „ ' 

Bquad I — , , 



2r4 



r4 



(9) 
(10) 

(11) 



respectively. We assumed a Bo nnor-E bert sphere (EES) den- 
sity profile (Bonnor [19561 Ebert |1955 ) with a central density of 
(^5) 2 X 10^^ m"^, a radius of the inner flat region of 500 AU, and 



an outer radius of 0.1 pc. We adopted a wavelength of 850 ^m, 
a constant dust temperature equal to 10 K for a source at the 
distance of 140 pc, and a resolution of 0'.'5. Images have then 
been smoothed to a beam of 10'' to get sharper contours. The 
magnetic field configurations are shown for four diff'erent incli- 
nations with respect to the plane of the sky, from zero (face-on 
source) ton/ 2 radians (edge-on source). 

Figure [T] refers to a purely toroidal magnetic field. As ex- 
pected, for a face-on source field lines are concentric while, ap- 
www.ita.uni-heidelberg.de/~dullemond/software/radmc-3qbroaching to the edge-on configuration, they become parallel in 



where i// is the angle between the north direction in the plane of 
the sky and the component of the magnetic field (B) in this plane 
and y is the angle between the direction of B and the plane of 
the sky. E and are given by 



f 

Jo 



BAiTdust) e dr^ 



(7) 
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the plane of the sky and can be pictured as concentric rings per- 
pendicular to the vertical axis of symmetry of the source. It can 
be easily understood that the depolarisation in the outer parts is 
due to the bending of the field lines which become perpendicular 
to the line of sight. 

Figures |2] and [3] show a radial configuration of the polarisa- 
tion vectors for a face-on source as well as the expected dipole 
and quadrupole field shape arising as the inclination increases. 




0.06 



0.04 



0.02 
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Fig. 1: Polarisation maps for a purely toroidal magnetic field 
and a BES density distribution (model radius, rmodei = 2.25 x 
10"^ AU). Black vectors represent the direction of the magnetic 
field, whose length is proportional to the fractional polarisation, 
obtained by flipping by 90 degrees the dust polarisation vectors. 
White contours show 50 and 90 per cent of the 850-yum dust 
emission (Stokes /) peak and the coloured map depicts the po- 
larised flux intensity (7^). The scale bar to the right gives Ip in 
Jy/(10'' beam). The inclination of the models with respect to the 
plane of the sky is given in the upper left corner of each subplot. 
The synthesised beam is shown in the lower left corner. 



4.2. Analytical models 

In order to test the correctness of our code, we reproduced the 
model described in Li & Shu (1996 ), hereafter LS96, to com- 
pare our results with those obtained by GGW05. This detailed 
comparison of the LS96 model between the results of GGW05 
and our code is worthwhile to have a further validation of the re- 
liability of the code. The model of LS96 is a semi-analytical, 
self-similar description of an axi-symmetric isothermal, self- 
gravitating protostellar core. It reproduces the hourglass mag- 
netic field configuration suggested by observations (Girart et 
al. 2006). The governing equations hinge on a dimension-less 
parameter, Hq, which characterises the fraction of support pro- 
vided against self-gravity by poloidal magnetic fields relative to 
that of gas pressure (see LS96 for more details). In particular. 
Fig. |4] shows the results for a core with Hq = 1.25, observed at 
850 yum. We adopted a distance of 140 pc, Tdust = 10 K and an 
angular resolution of O'.'l. Images have been smoothed to V for 
plotting. 




100 50 -50 -100 100 50 -50 -100 
Aa [arcsec] 

Fig. 2: Polarisation maps for a dipole magnetic field and a BES 
density distribution (rmodei = 2.25 x 10^ AU). See Fig. [l] for 
further information. 
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Fig. 3: Polarisation maps for a quadrupole magnetic field and 
BES density distribution (rmodei = 2.25 x 10"^ AU). See Fig.[l]for 
further information. 



Our computation supports the conclusions of GGW05 con- 
cerning the depolarisation at intermediate inclinations. In the 
central regions the bending of the field lines is stronger and the 
integration along line of sights close to the centre leads to a lower 
polarisation degree. We compared the minimum and the maxi- 
mum degree of polarisation as a function of the inclination of the 
core with values computed by GGW05. For a direct comparison, 
we smoothed our model results to 12'' which is the beam used 
in GGW05 and we obtain a substantial agreement as shown in 
Fig. in 

We examined the dependence of the polarisation degree 
on wavelength, temperature, and Hq, finding analogies with 
GGW05. The two top panels of Fig. [6] show the dependence of 
the polarisation degree on the wavelength. Due to the fact that 
the emission for lower wavelengths is rather flat and less concen- 
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Aa [arcsec] 

Fig. 4: Polarisation maps for the //q = 1.25 singular isothermal 
toroid according to the LS96 model (rmodei = 2.8 x 10^ AU). 
White contours show 10, 30, 50, 70, and 90 per cent of the 850- 
yum dust emission (Stokes /) peak and the scale bar to the right 
gives the polarised flux intensity in Jy/(r' beam). See Fig.[T]for 
further information. 
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20 40 60 80 
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Fig. 5: Maximum (solid lines) and mininum (dotted lines) per- 
centage of polarisation degree as a function of the inclination of 
the core with respect to the plane of the sky. DustPol model, 
solid squares; GGW05 model, empty squares. 



trated (Gon^alves et al. l2004K the region with low and uniform 
values of p, the so-called polarisation hole, has a smaller extent 
in this Figure at 450 jim with respect to 850 jim. It follows that, if 
the polarisation hole has a fixed physical radius, the correspond- 
ing value of ///max will be larger at 450 jim than at 850 jim. We 
also probed the dependence of the polarisation degree on the dust 
temperature as shown by the comparison of the two left panels 
of Fig. [6] The lower left panel shows the result for the assump- 
tion of an external heating with r^ust going from ~15 to ~8 K 
inside a radius of 0.1 pc (Gon9alves et al. [2Q04i) . As in GGW05, 



we conclude that a non-isothermal distribution of the dust tem- 
perature may help the observation of the depolarisation eff'ect. 
In fact, the magnetic field component in the plane of the sky is 
generally larger far from the central thermal dust emission peak. 
Thus, when in presence of an interstellar radiation field, the con- 
tribution of the external layers of the prestellar core to the Stokes 
parameters becomes substantial. Finally, we examined the de- 
pendence of the polarisation degree on //q, which is altering the 
magnetic field support. The lower right panel of Fig. [6] shows 
the distribution of the polarisation degree for //q = 0.125, corre- 
sponding to a mild pinching of the magnetic field lines. As found 
in GGW05, it is interesting to notice that even a low value of Hq 
leads to modifications of the polarisation degree for increasing 
intensities, measured as the diff'erence between the maximum 
(/'max) and the minimum (pmin) degree of polarisation observed. 
For Hq = 0.125, we found Ap = /^max - Pmin = 4.5% while, 
for comparison, Ap = 7.6% for Hq = 1.25, that is for a more 
pinched field, assuming in both cases v = 850 yum, r^ust = 10 K, 
and Qfmax = 15%. 




0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 



Fig. 6: Degree of polarisation (p) as a function of the intensity 
normalised to its peak value (///max) for the LS96 toroid with 
an inclination of 50 degrees with respect to the plane of the sky 
for diff'erent wavelengths, temperatures, and values of Hq. Each 
point in these diagrams represents a grid point from the model 
and multiple tracks are an artefact due to the grid sampling. The 
red dashed lines show the maximum and minimum values of the 
polarisation degree found by GGW05, see their Fig. 8. Regions 
on the right of the vertical blue solid lines correspond to the 
extent of the polarisation holes. 

The magnetic field geometry and temperature distribution 
are important components for the description of the observed 
values of the polarisation degree, but they cannot account for the 
whole set of observational constraints. For example, it is known 
that the maximum degree of polarisation, Ofmax, is strongly de- 
pendent on the cutoff' in the grain size distribution (Cho & 
Lazarian 20051) and depolarisation may be the result of alter- 
ations in the grain structure at higher densities which may reduce 
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the dust grain alignment efficiency (Lazarian & Hoang [20Q7]) . In 
order to test the dependence of the polarisation on grain prop- 
erties, we assumed a step function for Ofmax, which drops from 
15% to zero for molecular hydrogen number densities higher 
than 10^, 10^, 10^, and 10^ cm"^. Results are shown in Fig.[7] As 
expected, the polarisation degree decreases with intensity more 
rapidly as the limit on ^(H2) decreases. Despite the assumption 
of the step function is purely arbitrary, it is possible to deduce 
that the case of no polarisation for ^(H2) > 10^ cm"^ has to be 
discarded since there is no observational evidence that p is zero 
for / > 0.2/max. Although the theory of grain alignment is still 
not completely developed, the extreme flexibility in the parame- 
terisation of Ofmax allows the user to easily explore the parameter 
space. 



mentum vectors are initially aligned. We extracted four snap- 
shots from the RAMSES run at four times (see caption of Fig. [8]), 
chosen in order to predict the evolution of the polarisation since 
magnetic field lines get tangled as the core rotates. The out- 
put of the RAMSES run at these time samples was used as in- 
put to DustPol using a distance of 140 pc, with a wavelength 
of 850 yum, Tdust - 10 K, and an angular resolution of 0^.^05, 
then smoothed at 0'.'5, which is a resolution typically achievable 
by ALMA. Note that these maps are not corrected for the in- 
strumental response of the interferometer. Initially, the magnetic 
field has hourglass shape, as time increases, the field configura- 
tion becomes purely toroidal. Figure [8] shows the time evolution 
for a face-on configuration, which yields an initially radial pat- 
tern changing to a more toroidal structure. 
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Fig. 7: Degree of polarisation {p) as a function of the intensity 
normalised to its peak value (I/Imax) for the LS96 toroid as in 
Fig. [6] This plot shows the decrease in p for high values of ///max 
due to the assumption of a step function for the maximum degree 



of polarisation. Labels (10^ to 10^ 
^(H2) above which ofmax = 0- 



cm represent the limit on 
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0.05 
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Fig. 8: Polarisation maps for a collapsing magnetised dense 
core (face-on view) obtained with the RAMSES code (rmodei = 
6.3 X 10^ AU). Labels refer to diff'erent times of the collapse: 
ti = 1.76 X 10^ yr, t2 = 1.86 X 10^ yr, ^3 = 1.94 X 10^ yr, and 
t4 = 2.96 X 10"^ yr. White contours show 30, 50, 70, and 90 per 
cent of the 850-yum dust emission (Stokes /) peak and the scale 
bar to the right gives Ip in Jy/(0'.'5 beam). See Fig. [T] for further 
information. 



4.3. Numerical models 

The novelty of DustPol is that it can easily manage inputs from 
numerical simulations set out on regular or irregular grids. In 
this paper we give a couple of examples of 3D numerical cal- 
culations of a coll apsing magnetised den se cor e (Hennebelle 
& Fromang l20Q8l Hennebelle & Ciardi [20091) which make 
use of the adaptive mesh refinement (AMR) code RAMSES 
(Teyssier 2002, Fromang et al. 2006). We describe only quali- 
tatively the results postponing quantitative conclusions for a fol- 
lowing paper (Padovani et al. in prep.). 

The grid of a numerical simulation is sampled as if it were an 
analytic function. Each point of the grid built by LIME takes the 
value of density, temperature, and magnetic field of the closest 
model point. Before this step, the output of a numerical simula- 
tion has to be converted in the format readable by LIME. To this 
end, a simple wrapper is distributed together with DustPol. 

RAMSES was run for a collapsing rotating core with mass- 
to-flux ratio equal to 5. The magnetic field and the angular mo- 



Finally, we tested DustPol using a large scale RAMSES 
simulation (Dib et al. 2010) relative to a prestellar core 
with a mass of 1.62 M©, which has been obtained from the 
"Molecular cloud evolution with decaying turbulence' 'project 
of the StarFormat databas^ The upper panel of Fig. |9[ shows 
the column density distribution of molecular hydrogen, A/^(H2), 
of this core. DustPol has been used to get the polarisation pat- 
tern using a distance of 140 pc, r^ust = 10 K, and a wavelength 
of 850 yum. The resulting image is shown in the lower panel of 
Fig.|9] It has been smoothed from V to 14V5 which is the beam 
of SCUBA, accounting for the error beams for an accurate de- 
scription of the beam profile (Hogerheijde & Sandell 2000 ). 

Under the optically thin assumption, the agreement between 
the RAMSES H2 column density distribution and the DustPol 
Stokes / map is remarkable. In fact, after normalising N(}l2) and 
Stokes / with respect to their maximum values, we computed the 
ratio between N(ll2) and Stokes / finding a mean value of 0.95 ± 

^ http : //starf ormat . obspm . f r/starf ormat/pro j ects 
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0.12. Additionally, it is worth noting the displacement between 
the polarised emission and the dust continuum emission near the 
offset (0, 0) and, conversely, the good overlap near (-100, 150). 
This can be easily understood by looking at the local direction of 
the magnetic field into the RAMSES data cube. We found that 
the magnetic field points towards the line of sight around the 
centre of the map, thus explaining the decrease in polarisation 
degree. On the other hand, the direction of the magnetic field 
mainly lays on the plane of the sky near the north-west clump, 
where the polarisation is stronger. 

Working in tandem, DustPol and StarFormat (to be publicly 
available) can be used to interpret large-scale polarisation maps 
of molecular clouds, investigating on magnetic field morpholo- 
gies. 




-100 



100 -100 

Aa [arcsec] 

Fig. 9: Upper panel, column density distribution of molecular 
hydrogen for a prestellar core extracted from the StarFormat 
database. Black and white contours show 30, 50, 70, and 95 per 
cent of the H2 column density peak and of the 850-yum dust emis- 
sion (Stokes /) peak, respectively. The scale bar to the right gives 
A^(H2) in logarithmic scale. Lower panel: polarisation map for 
the core depicted in the upper panel (rmodei = 2.8 x 10"^ AU). 
The scale bar to the right gives Ip in Jy/(14V5 beam) and it refers 
to the coloured map which has been debiased for values below 
3(T. Magenta contours show 30, 50, 70, and 95 per cent of the 
850-yum dust emission (Stokes /) peak; white vectors show the 
direction of the magnetic field. The synthesised beam is shown 
in the lower left comer. 



4.4. Comparison with observations: NGC 1333 IRAS 4A 

The ultimate goal of the DustPol project is to generate simula- 
tions capable of being compared to the real data. Its main advan- 
tage is that it provides quick convergence because of the adaptive 
gridding and it easily accommodates outputs from other codes 
as shown in Sect. |4.2| and [4.3[ Recently, FGGll have devel- 
oped a method to compare the observational polarisation maps 
of the low-mass Class object NGC 1333 IRAS 4A to a set 
of magnetohydrodynamic (MHD) simulations, which allows to 
find the best set of models that fit the magnetic field properties 
of NGC 1333 IRAS 4A. They generate a synthetic source based 
on models and simulations and, using the same equation set de- 
scribed in this work, they provide synthetic Stokes maps. These 
synthetic maps are convolved by the interferometric response 
in order to obtain the spatially filtered maps comparable to the 
real data. To assess the quality of the DustPol module, we re- 
peated this scheme replacing the radiative transfer tool in FGGl 1 
with DustPol. We used the Allen et al. (2003 ) model datacube 
with the best fitting parameters found by FGGl 1, namely a mass 
of 1.6 Mo, a mean volume density within the angular size of 
100 AU of 3 X 10^ cm"^, an inclination of 45 degrees with re- 
spect to the plane of the sky, Hq = 0.125 (see Sect. |4.2| ), and the 
dust temperature profile derived by Maret et al. ( 2002 ). We ran 
DustPol with a resolution of 0'.'16, assuming a source distance 
of 300 pc, and a wavelength of 880 jim. Then, we applied the 
same convolution tool developed first by Gongalves et al. (20081 
and improved by FGGll to the output map of DustPol for a 
final synthesised beam of IV 24 x IV 12. Figure 10 shows our re- 
sults, to be compared with the third row of Fig. 8 in FGGll. 
The observations and the model are in good agreement and the 
general morphologies of the total emission (Stokes /) and the po- 
larised emission are consistent. In fact, even if the polarisation 
vectors in the map generated by DustPol cover a larger area, the 
diff'erence map between the Stokes Q (and U) maps created by 
the two radiative transfer tools show no features above the 3cr 
level (cr = 2.5 mJy beam"^). Therefore, the diff'erences are due 
to the statistical noise at the low RMS regions. Finally, the dif- 
ferences between the position angles evaluated with the radiative 
transfer tool in FGGll and DustPol have a standard deviation 
of 0.48 degrees, confirming the reliability of DustPol. 



5. Conclusions 

We have presented the DustPol module of the Adaptable 
Radiative Transfer Innovations for Submillimetre Telescopes 
(ARTIST) package, a tool for modelling the polarisation of ther- 
mal dust emission self-consistently. The module is an exten- 
sion to the radiative transfer code LIME making it possible to 
model complex three-dimensional structures. The synergy be- 
tween DustPol and codes such as RAMSES and RADMC-3D 
will help to enlighten the interpretation of the role of magnetic 
fields in star-forming regions. 

DustPol can easily handle analytical as well as pre-gridded 
models and outputs from other codes, accounting for the depen- 
dence on dust temperature and, for the first time to our knowl- 
edge, for variations of the polarisation degree. A number of ex- 
amples (Sect. [4]) demonstrate the capabilities and the accuracy 
of DustPol to interpret existing submillimetre data as well as 
new data from ALMA and the Planck satellite. The output of 
DustPol is a FITS file whose header is consistent with the for- 
mat of the input models for the sim_observe-sim_analyze 
tasks of the CAS A package. This makes it straightforward, e.g.. 
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Fig. 10: Top panel: dust polarisation SMA maps toward 
NGC 1333 IRAS 4A from FGGll. Middle panel: synthetic 
dust polarisation map generated with DustPol, following the 
FGGl 1 convolution scheme, using the best fitting parameters of 
the Allen et al. (2003) model found by FGGll (see Sect. [44 
for more details). Each panel shows the intensity (contours), 
polarised intensity (pixel map), and direction of magnetic field 
(blue segments). Contours depict emission levels from 6(T up to 
the maximum value in steps of 6cr, where a = 0.02 Jy beam"^ 
The middle panel shows the beam and the angular and spatial 
scale. Bottom panel: histogram of the diff'erence between the po- 
sition angles computed with the radiative transfer tool in FGGl 1 
and DustPol. 



to make prediction about the polarised flux for a given configu- 
ration of the ALMA antennas. 

The code is publicly available as long as appropriate refer- 
ence is made to this paper. It is accessible through the complete 
ARTIST package, taking advantage of the Python-based graphic 
user interface that includes a model library or through a stand- 
alone version. 
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